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Abstract 

We continue a previous work on the comparison between the post-Newtonian (PN) approximation 
and the gravitational self-force (SF) analysis of circular orbits in a Schwarzschild background. We 
show that the numerical SF data contain physical information corresponding to extremely high PN 
approximations. We find that knowing analytically determined appropriate PN parameters helps 
tremendously in allowing the numerical data to be used to obtain higher order PN coefficients. 
Using standard PN theory we compute analytically the leading 4PN and the next-to-leading 5PN 
logarithmic terms in the conservative part of the dynamics of a compact binary system. The 
numerical perturbative SF results support well the analytic PN calculations through first order in 
the mass ratio, and are used to accurately measure the 4PN and 5PN non-logarithmic coefficients 
in a particular gauge invariant observable. Furthermore we are able to give estimates of higher 
order contributions up to the 7PN level. We also confirm with high precision the value of the 3PN 
coefficient. This interplay between PN and SF efforts is important for the synthesis of template 
waveforms of extreme mass ratio inspirals to be analysed by the space-based gravitational wave 
instrument LISA. Our work will also have an impact on efforts that combine numerical results in 
a quantitative analytical framework so as to generate complete inspiral waveforms for the ground- 
based detection of gravitational waves by instruments such as LIGO and Virgo. 
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I. MOTIVATION AND SUMMARY 



This paper is the follow up of previous work [T] (hereafter Paper I) where we demonstrated 
a very good agreement between the analytical post-Newtonian (PN) approximation and the 
numerical gravitational self-force (SF) for circular orbits in the perturbed Schwarzschild 
geometry. The first step had been taken by Detweiler [2] who showed agreement at 2PN order 
using the existing PN metric [3]. 1 Motivated by this result we pushed the PN calculation in 
Paper I up to the 3PN level. This is particularly interesting because the 3PN approximation 
necessitates an extensive use of dimensional regularization to treat the divergent self-field 
of point particles. The successful comparison reported in Paper I confirmed the soundness 
of both the traditional PN expansion (see e.g. |4]) and the perturbative SF analysis [5]- 
IH] in describing the dynamics of compact binary systems — notably, regarding subtleties 
associated with the self-field regularizations in use in both approaches. This comparison 
dealt with the conservative part of the dynamics, but previous comparisons between the PN 
and the SF had checked dissipative effects [TOHTT] . 

In Paper I we also showed that the quality of the numerical SF data is such that sub- 
stantial physical information remains far beyond 3PN order, i.e. is contained within the 
numerically derived residuals obtained after subtracting the known 3PN terms from the 
data (see Fig. 3 of Paper I). In the present paper we explore further the higher-order PN 
nature of the numerical data. We point out that knowing analytically determined appropri- 
ate PN parameters helps tremendously in allowing our numerical data to be used to obtain 
higher order PN terms. In particular, we show that prior analytic information from PN the- 
ory regarding the presence of logarithmic terms in the PN expansion is crucial for efficiently 
extracting from the SF data the numerical values of higher order PN coefficients. 

The occurence of logarithmic terms in the PN expansion has been investigated in many 
previous works [T8H25] . Notably Anderson et al. [21] found that the dominant logarithm 
arises at the 4PN order, and Blanchet & Damour |25j (see also [26J) showed that this 
logarithm is associated with gravitational wave tails modifying the usual 2.5PN radiation- 
reaction damping at the 1.5PN relative order. Furthermore the general structure of the 
PN expansion is known [21]: it is of the type Y2( v /c) [ln(u /c)] q , where k and q are positive 
integers, involving only powers of logarithms; more exotic terms such as [ln(ln(t> /c))] q cannot 
arise. In the present paper we shall determine the leading 4PN logarithm and the next-to- 
leading 5PN logarithm in the conservative part of the dynamics of a compact binary system. 

Consider two compact objects with masses mi and (without spins) moving on ex- 
actly circular orbits. The dissipative effects associated with gravitational wave emission are 
neglected, which is formalized by assuming the existence of a helical Killing vector field 
K a (x), being null on the light cylinder associated with the circular motion, time-like inside 
the light cylinder (for instance at the particle's location) and space-like outside (including a 
neighborhood of spatial infinity). Then we consider a particular gauge invariant observable 
quantity [2] defined as the constant of proportionality between the four-velocity of one of 
the masses, say mi, and the helical Killing vector evaluated at the location of that particle, 
i.e. K? = K a ( yi ), 

u1=ulK«. (1.1) 



As usual the nPN order refers to terms equivalent to (v/c) 2n beyond Newtonian theory, where v is a 
typical internal velocity of the material system and c is the speed of light. 
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The quantity u± represents the redshift of light rays emitted from the particle and received 
on the helical symmetry axis perpendicular to the orbital plane [2]; we shall sometimes refer 
to it as the redshift observable. Adopting a coordinate system in which the helical Killing 
vector field reads K a d a = d t + Vtd^, where Q denotes the orbital frequency of the circular 
motion, we find that the redshift observable reduces to the t component u\ = dt/dri of the 
particle's four-velocity, namely 

^ = u\=(-g aP ( yi )^j V2 . (1.2) 

Here = dyf/dt = (c, v\) is the ordinary coordinate velocity used in PN calculations, 
and g a p(yi) denotes the metric being evaluated at the particle's location by means of an 
appropriate self-field regularization, i.e. mode-sum regularization in the SF approach, and 
dimensional regularization in the PN context. 

The point is that uj can be computed as a function of the orbital frequency Q in both the 
PN approach for any mass ratio, and in the perturbative SF framework when the mass mi is 
much smaller than m-i- Summarizing the analytical 3PN result of Paper I and present com- 
putation of the 4PN and 5PN logarithmic terms in Sees. [TI}fV| we obtain the SF contribution 



to the redshift observable (1.2) as 



t4 = - y - 2y 2 - 5y 3 + ( ~ + ^vr 2 ) y 4 

( 64 \ / 956 \ 
+ ( a 4 - y lny) y 5 + f a 5 + — lnj/J y 6 + o(y 6 ) , (1.3) 

where y = (Gm 2 fi/c 3 ) 2//3 is a PN parameter associated with the lighter mass m 1; and 
and «5 denote some purely numerical coefficients left out in the PN calculation. However, 
having obtained theoretical predictions for the 4PN and 5PN logarithmic terms, we are able 
to perform an efficient fit to the numerical SF data and to accurately measure the other non- 
logarithmic 4PN and 5PN coefficients. We find a 4 = -114.34747(5) and a 5 = -245.53(1) 
where the uncertainty in the last digit is in parenthesis. Furthermore we can also measure the 
6PN coefficients and /3q (such that + (5q In y is the factor of y 7 ), and give an estimate of 
the total contribution of the 7PN coeffic ient (i ncluding both logarithmic and non-logarithmic 



V] and Fig. [l] in Sec. f 
reement with the SF d 



IVID The 3PN coefficient a 3 = + %tt 2 is also 



terms); see Table 

found to be in agreement with the SF data with high precision. 

The non-logarithmic coefficients a^, a$, ■ ■ ■ would be extremely difficult to obtain with 
standard PN methods. Their computation would require in particular having a consistent 
self-field regularization scheme; for instance it is not guaranteed that dimensional regular- 
ization which has been so successful at 3PN order could be applied with equal success at 
much higher orders. Nevertheless these coefficients are obtained here for the first time with 
reasonable precision up to the impressive 7PN order. This emphasizes the powerfulness of 
the perturbative SF approach and its ability to describe the strong field regime of compact 
binary systems, which is inaccessible to the PN method. Of course, the limitation of the SF 
approach is the small mass-ratio limit; in this respect it is taken over by the PN method. 

The analytical and numerical results obtained in this paper up to 7PN order could be 
used for the synthesis and calibration of template waveforms of extreme mass ratio inspirals 



Inspired by our earlier work [T] , the easy calculation of the 4PN logarithm has already been given in [57] . 
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to be observed by the space-based gravitational wave detector LISA. They are also relevant 
to analyses that combine numerical computations in a quantitative analytical framework for 
the generation of inspiral waveforms for the ground-based LIGO and Virgo detectors. 

The remainder of this paper is organized as follows: In Sec. [IT] we perform a detailed 
analysis of the occurence of logarithmic terms in the near-zone expansion of an isolated 
source. This general discussion is followed in Sec. Ill by the explicit computation of the 
leading order 4PN and next-to-leading order 5PN logarithmic terms in the near-zone metric 
of an arbitrary post-Newtonian source, and then of a compact binary system. We proceed 
in Sec. IV with the computation of these terms in the acceleration of the compact binary, 
as well as in the binary's conserved energy, and consider the restriction to circular orbits. 
This allows us to derive intermediate results necessary for the computation of the 4PN and 
5PN logarithmic terms in the redshift observable (1.2) for circular orbits; this is detailed in 
Sec. |Vj Finally, Sec. VI is devoted to a high-order PN fit of our numerical data for the SF 
effect on the redshift variable. The Appendix provides general formulas for the computation 
of logarithmic terms in PN theory. 



II. GENERAL STRUCTURE OF LOGARITHMIC TERMS 



A. Near-zone expansion of the exterior metric 



In this Section we study in a general way the PN orders at which logarithmic terms 
occur in the near-zone expansion of the metric of an isolated source. Our main tool will be 
the multipolar-post-Minkowskian (MPM) analysis of the vacuum field outside the compact 
support of the source [24Tl2"6| |28| [29] . The starting point is the general solution of the 
linearized vacuum Einstein field equations in harmonic coordinates, which takes the form 
of a multipolar expansion parametrized by mass-type Ml and current-type Sl multipole 
moments 1 30 1 3 
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3 Our notation is as follows: L — i\ - ■ ■ %i_ denotes a multi-index composed of £ multipolar spatial indices 
ii,--- ,ie (ranging from 1 to 3); &l — d% x ■ ■ ■ di e is the product of i partial derivatives di = d/dx l ; 
x l = x i t " " %i e is the product of £ spatial positions xf, similarly = n il ■ ■ ■ rii e is the product of I unit 
vectors rii = Xi/r; the symmctric-trace-free (STF) projection is indicated with a hat, i.e. xl = STF[xl], 
fiL = STF[ni], &l = STF[<9l], or sometimes using brackets surrounding the indices, i.e. x^l) = xl- In 
the case of summed- up (dummy) multi-indices L, we do not write the £ summations from 1 to 3 over 
their indices. The totally antisymmetric Levi-Civita symbol is denoted £ijk', symmetrization over indices 
is denoted (ij) = + ji); time-derivatives of the moments are indicated by superscripts (n). 
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The multipole moments Ml and Sl are symmetric and trace-free (STF) with respect to all 
their indices and depend on the retarded time u = t — r/c in harmonic coordinates. They 
describe a general isolated source and are unconstrained except that the mass monopole M 
and current dipole Si are constant, and the mass dipole Mj is varying linearly with time. 

Starting from hi we define a full non-linear MPM series for the "gothic" metric deviation 
h af5 = yj—g g a ^ — i] a P (where g a ^ and g denote the inverse and determinant of the usual 
covariant metric respectively, and where rj a ^ is the Minkowski metric) as 

+00 

h"? = G n hf , (2.2) 

n=l 

where the Newton constant G serves at labelling the successive post-Minkowskian orders. 
Plugging this series into the (vacuum) Einstein field equations in harmonic coordinates we 
find at each order d^h"^ = 0, together with 

Uhf = Nf , (2.3) 

where □ = rj^d^d^ is the flat d'Alembertian operator, and where N n denotes the n-th 
non-linear gravitational source term depending on previous iterations hi, ■ ■ ■ , h n _i. An 



explicit "algorithm" has been proposed in [21] for solving (2.3) and the condition of harmonic 
coordinates at any post-Minkowskian order n. 



We are interested in the expansion of the solution of (2.3) in the near-zone (NZ), i.e. 
formally when r — » (but still outside the compact supported source). The general structure 
of that expansion is known [26]. For the source term we have (the NZ expansion being 
indicated with an overbar) 



We see that besides the normal powers of r we have also powers of logarithms of r; p is 
an integer (p G Z) bounded from below by some po depending on E n , and q is a positive 
integer (q G N). We pose A = 27rc/f2, with Q a typical frequency scale in the source 
to be identified later with the orbital frequency of the binary's circular orbit. We denote 
by E n = {Mlu Ml 2 , ■ ■ ■ 1 £ait +1U SaLn-i} a se t °f n multipole moments, with the current 
moments endowed with their natural Levi-Civita symbol. We pose £ t = li for mass moments 
and li = li + 1 for current moments, so that Yl7=i £i ^ s ^ e total number of indices carried 
by the moments of the set E n . On the other hand I is the number of indices carried by the 



STF multipolar factor fit,. The multipole functions in (2.4) admit the general structure 



F&V) = / dMl '"/ du n )Cf Li ... Ln (t,ui,--- ,u n )M^\ui)---e aien+lle XLl-iM, (2.5) 

where the kernel K, has an index structure made only of Kronecker symbols and is only a 
function of time variables: the current time t, the n integration arguments u, (satisfying 
Ui ^ t) and the period P = X/c of the source. Then with this convention we see that the 



powers of 1/c in (2.4) are set by dimensionality. A useful lemma [26] is the fact that the 



multipolar order I is necessarily constrained by the following two inequalities: 

n n 
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Here s is the number of spatial indices among a and /?, i.e. the "spin" given by s 
according to a (3 = 00,0i,ij. 



0,1,2 



The lemma (2.6 ) will serve at controlling the PN order of "branches" of logarithmic terms 
arising in the MPM iteration of the external field. Already we know [21] that the powers 
of the logarithms are limited to q ^ n — 2 in the source term N n . After integration of the 
source term N n we shall find the corresponding solution h n which will admit the same type 
of NZ expansion as its source. However the maximal power of the logarithms in the solution 
will be increased by one unit with respect to the source and is thus limited by n — 1, i.e. 
q ^ n — 1 in h n . For instance this means that logarithms squared cannot arise before the 
cubic non-linear order n — 3. 

To control the occurence of logarithms in the near- zone it will be sufficient to integrate 
the source (2.4) by means of the integral of the "instantaneous" potentials defined by formal 
PN iteration of the inverse Laplace operator A -1 , say = A -1 + c~ 2 <9 2 A~ 2 + ■ ■ • . This 
is because any homogeneous solution to be added to that particular solution will have 
the structure of a free multipolar wave (retarded or advanced) whose near-zone expansion 
cannot contain any logarithms. However, when acting on a multipolar expanded source 
term, valid only in the exterior of the matter source and becoming singular in the formal 
limit r — > 0, we must multiply the source term by a regulator (r/A) B , where B is a complex 
number and A = cP is the length scale associated with the orbital motion. After applying 
the instantaneous propagator we take the finite part (FP) of the Laurent expansion when 
B — > 0. Thus the solution reads as 



-a/3 



+oo 



fp y 

B=0 



k=0 



d_ 
cdi 



2 k 



-k-l 



\) n 



(2.7) 



Later, in (2.16) below, we shall denote by I 1 the particular "instantaneous" regularized 
propagator appearing in (2.7). The term H n denotes the NZ expansion of an homogeneous 



solution of the d'Alembert equation. In the general case this solution will be a mixture of 
purely retarded and advanced multipolar waves, say of the type ^0l{.Rl(£ — r/c)/r} and 
Y^dL{A L {t + r/c)/r}, but the point is that the NZ expansion of H n when r — > clearly 
does not contain any logarithms. So in order to control the logarithms we can ignore the 
homogeneous piece H n . 

As argued in [26] the use of the latter "instantaneous" propagator, say X -1 , corresponds 
to keeping only the conservative part of the dynamics, i.e. neglecting the dissipative part 
associated with gravitational radiation-reaction. Below we shall implement the restriction to 
the conservative case by looking at circular orbits with helical Killing symmetry. We expect 



that a solution admitting this symmetry should be given by (2.7) where the homogeneous 
part H n is of the symmetric type ^^{[Sz^t — r/c) + S^t + r/c)]/r}. In this "symmet- 
ric" situation, where the radiation-reaction is neglected, the solution should depend on the 



length scale A appearing in the first term of (2.7). Indeed this length scale is introduced in 



the problem by our assumption of having the helical Killing symmetry with Killing vector 
K a d a = dt + Vtdy where Vt = 2nc/ A. 



B. Near-zone versus far-zone logarithms 



Inserting the general form of the source term (2.4) into (2.7), and ignoring from now 
on the homogeneous term H n which does not contain logarithms, we obtain (dropping the 
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space-time indices a(3 for clarity) 
1 



h, 
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We can explicitly integrate the iterated Poisson integral and find 



A 
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with S-dependent coefficients 



(2.* 



(2.9) 



' P ' fc(5) El (B+p + 2 + 2i-£)(B+p + 3 + 2i + 



(2.10) 



We shall now control the occurence of a pole oc 1/B in the latter expression which, after 
taking the finite part in (2.8), will generate a logarithm lnr. Actually, since we have to 
differentiate q times with respect to B, the pole in a^ tP ^{B) (which is necessarily a simple 
pole) will yield multiple poles oc 1/B m , and we shall finally end up with powers of logarithms 
(lnr) m , where here m ^ q+ 1 — hence the increase by one of the powers of logarithms from 
the source to the solution, as discussed previously. 

Inspection of Eq. (2.10) readily shows that there are two types of poles. First we have 
the poles for which p + 2 = £ — 2i. These will be qualified as "near-zone poles", and the 
structure of the solution for these poles reads 



NZ pole 
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(2.11) 



where j = k — i ^ and the functions Gi,j m (t) have a structure similar to (2.5). Note that 
(2.11) is perfectly regular when r — > [at least when £ + j ^ 1] and will therefore be valid 
(after matching) inside the matter source. On the other hand the "far- zone poles" for which 
p + 2 = — £ — 1 — 2i have the structure 



FZ pole 
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t,i>o 
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r\ 2j 
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(2.12) 



These poles become singular when r — > 0. We shall argue later that the associated logarithms 
do not contribute to the PN expansion of quantities we compute in this paper (like the 
redshift observable or the conserved energy of a compact binary system). 

We can now easily control the PN order of these poles. Taking into account all the powers 
of 1/c and the fact that j ^ 0, we obtain 



(h n ) 



NZ pole 



E° 



c 3 ™+Er=i^ 



(2.13) 



Next, the inequality in the left of the lemma (|2.6|) provides a uniform bound of the PN order 

1 



of each of the terms in (2.13), leading to 

(hn) 



NZ pole 
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(2.14) 
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This means that the NZ poles in the n-th non-linear metric are produced at least at the 
^y^PN level; note that the power of 1/c in the ij components of the perturbation h n , such 
that s = 2, gives immediately the PN order. Similarly we find, using now the inequality on 



the right of (2.6), that the FZ poles are produced at the level 



(K) 



FZ pole 



o 



■,3n—l—s 



(2.15) 



corresponding to (at least) the ^^y^PN order. Notice that the far-zone poles come earlier 
than the near-zone ones in the PN iteration. 

We use these general results to control the occurence of (powers of) logarithms in the PN 
expansion. First be careful that our findings do not mean that all the logarithms at some 



n-th non-linear order will have the PN orders (2.14) and (2.15); it states that whenever new 



logarithms appear they are necessarily produced at least at these PN levels. However, once 
a "new" logarithm has been produced in h n , it will contribute in the source term N n+ i of 
the next iteration, and therefore will also appear in the corresponding solution h n+ \ where 
it needs not to be associated with a pole occuring at that order. In fact we expect that the 
vast majority of logarithms only come from the iteration of original logarithms seeded by 



poles. Such "iterated" logarithms will escape the rules (2.14) and (2.15). 



Given a logarithm at order n coming from a NZ pole and being thus at least of order 
^Y^PN, we can check that it will generate iterated logarithms at any subsequent non-linear 
order n + p, with p ^ 1, and that those will be at least of order 3n 2 2p PN. We can therefore 
always bound the PN order of the complete family of iterated NZ logarithms by the order 
3j y^PN of the "seed" logarithm. 4 The same reasoning applies for the PN orders of the 
iterated FZ logarithms which are bounded from below by the ^y^PN order of the seed. 

When n = 2 we find from (2.14) that there is a family of NZ logarithms starting at the 



4PN order. We know that the 4PN logarithmic term is associated with gravitational wave 
tails; it has been computed for general matter sources in [25]. Conjointly with this 4PN 
logarithm there will be also logarithms at 5PN and higher orders, all of them at quadratic 
order n = 2, and all these quadratic logarithms will have to be iterated at the next cubic 
order n — 3, and so on. As we discussed this defines a complete family of NZ logarithms, 
and this family will be sufficient to control all the NZ logarithms at 4PN and 5PN orders. 
Indeed, we expect that at cubic order n = 3 a new family of NZ logarithms will appear, but 
according to the result (2.14) this new family will be of order 5.5PN at least. In particular 
this reasonning shows that the dominant NZ logarithm squared [ln(r/A)] 2 is at least 5.5PN 
order. Such 5.5PN logarithm would be time odd in a time reversal and belongs to the 
dissipative radiation-reaction part of the dynamics so we shall ignore it. Similarly the next 
family coming at the quartic approximation n = 4 will be at least 7PN — thus the dominant 
[ln(r/A)] 3 is expected to appear at least at 7PN order. 

We shall now argue that only the family of NZ logarithms starting at the 4PN order needs 
to be considered for the present computation, because quite generally the FZ logarithms 
cannot contribute to the conserved part of the dynamics of a compact binary system. 



When p = 1 we get the same PN order as the seed logarithm because according to ( 2.14 1 the ij component 
of the metric perturbation is of order 1 /c 3n+2 , and hence generates at the next iteration a term of order 
l/c 3rl+4 in the 00 component of the metric perturbation h"+i (via the non-linear source term h l Jdidjhf), 
which is still of 3 "+ 2 PN order. We shall use later the trick that by gauging away the ij component of the 
metric perturbation we can greatly simplify the computation of the subsequent iteration. 
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C. Argument that far-zone logarithms give zero contribution 



(2.15) 



The FZ logarithms are generated by seeds whose PN order is controlled by the estimate 
First one can check that due to the particular structure of the quadratic metric 
n = 2 there is no FZ pole at the quadratic order [25J. The FZ logarithms come only at the 
cubic order n = 3 and from the estimate (2.15) we see that they arise dominantly at 3PN 
order, i.e. earlier than the NZ logarithms at 4PN order. The 3PN far-zone logarithms have 
been investigated in [25] and also in previous work [2T] . However we do not need to consider 
these and other FZ logarithms in the present calculation as the following argument shows. 
The NZ and FZ logarithms were investigated using the operator of the "instantaneous" 



potentials defined by [see Eq. (2.7) 
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(2.16) 



This propagator depends on the length scale A. Now our basic assumption is that in order 
to treat the conservative part of the dynamics, admitting the helical Killing vector K a d a = 
dt + Qdp in the two-body case, one should integrate the field equations with the propagator 
(2.16) in which we set A = 2itc/Vl. In this way the conservative dynamics will fundamentally 



depend on the scale A coming from the Killing symmetry and explicitly introduced through 



the propagator (2.16). 



By contrast, in a physical problem where we look for the complete dynamics including 
both conservative and dissipative (radiation-reaction) effects, there is no preferred scale such 
as A — indeed, nothing suggests that the dynamics should depend on some pre-defined scale 
A. In this case we integrate the field equations using the standard retarded integral, i.e. 



[Nn] 



FPZI 

B=0 47T 



d 3 



X 



l x 
T 



AL.(x',t 



(2.17) 



The non-linear source term N n is in unexpanded form since we integrate in all the exterior 



of the source and not only in the NZ as in (2.16). But, as in (2.16), we have introduced a 
regulator (|x'|/A) s and a finite part to cure the divergencies of the multipole expansion at the 



origin of the coordinates. Because of this regulator, the retarded integral (2.17) depends on 



the scale A which must therefore be cancelled by other terms in the physical metric. What 
happens is that the dependence on A coming from integrating the non-linearities using 



(2.17) is cancelled by a related dependence on A of the multipole moments of the source 



which parametrize the linear (retarded) approximation. The source multipole moments can 
be written as integrals over the pseudo stress-energy tensor of the matter and gravitational 
fields [28]. Because of the non-compactness of the gravitational field the integral extends 
up to infinity and involves a similar regulator (\x!\/X) B dealing with the boundary of the 
integral at infinity. The final independence of the physical metric on A can be checked by 
formally differentiating the general expression of the metric found in [28] . The cancellation 
of A has been explicitly verified up to the 3PN order in the case of compact binaries |31j . 
What is the difference between the physical situation and the "unphysical" one in which 



we would use the propagator (2.16) ? To compare the two situations we expand the retarded 



integral (2.17) in the near-zone. Recalling that the overbar refers to the NZ expansion, we 
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obtain 



□iW - 1- 1 F.] + E ^ " rh) ~J l(t + r/c) } ■ 

showing that the two solutions differ by an homogeneous solution of the wave equation 
which is of the anti-symmetric type (i.e. retarded minus advanced) and is therefore regular 
in the source. We know that the multipolar functions Tl(u) parametrizing this solution 
are associated with non-linear tails and their expressions can be found in [26, 29J. In the 



physical case, the homogeneous solution in (2.18) will remove the A-dependence located in 



the NZ logarithms appearing from the first term, and which have the symbolic NZ structure 
~ £z,ln(r/A). On the other hand the A-dependence in the FZ logarithms ~ du(l/r) ln(r/A), 
is removed by the retarded homogeneous solution we start with at the linear approximation. 
Now in the unphysical situation we shall want to subtract the anti-symmetric solution in 



(2.18) in order to use the instantaneous propagator X~ x . Therefore the scale A will no longer 
be cancelled from the near-zone logarithms ~ x^XivirjX) which will thus remain as they 
are. Suppose that they are evaluated at the location of a body in a two-body system, then 
the NZ logarithms become ~ y^ln(|yi|/A) where yi is the position of the body, and hence 
~ yi In(ri2/A) in the frame of the center of mass, where r^i is the two-body's separation. 
Using Kepler's law the logarithm becomes In(ri2/A) = ^ln7 where 7 = Gm/(ri2C 2 ) is a 
standard PN parameter, showing that the NZ logarithms do contribute to the final result. 

On the contrary the FZ logarithms ~ ^(l/r) ln(r/A) will not. Indeed the scale A therein 
will still be cancelled out by the linear retarded solution. 5 This means that in the application 
to binary systems the final FZ logarithms are scaled not by A but rather by the size r 12 of 
the orbit, and become some ~ di(l/r) ln(r/r 12 ). When considered at the location of one of 
the bodies we get ~ c?i(l/|yi|) ln(|yi|/ri 2 ) which clearly does not contribute in the center- 
of-mass frame. The latter reasoning is rather formal because the multipole expansion is 
valid only outside the source and it does not a priori make sense to apply it "at the location 
of one particle." However the reasoning may be better justified from a matching argument 
suggesting that the multipole expansion is valid "everywhere" , in a restricted sense of formal 
asymptotic series. 

Our conclusion is that we do not need to consider the FZ logarithms. From the previous 
investigation we see that it is sufficient to consider the family of iterated NZ logarithms 
generated at the quadratic order n = 2, and to compute the 4PN and 5PN logarithms 
within this family. We devote the next Section to this task. 



III. THE 4PN AND 5PN NEAR-ZONE LOGARITHMS 
A. External near-zone post-Newtonian metric 

Following [251 EE] we know that the dominant logarithms in the near-zone metric are 
coming from "tails" generated by quadratic coupling between the constant total mass M of 



5 The argument could be extended to an unphysical solution which would be truly symmetric in time, i.e. 
which would start with a symmetric (retarded plus advanced) linear approximation and integrate the 
non-linearities by means of the propagator I , 
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the system (i.e. the ADM mass) and the time varying multipole moments Ml or Sl- Let 



us define (n, u) as being the coefficient of the leading 1/r piece in the non-stationary or 



"dynamical" part (h* )dyn of the linearized metric given by (2.1), i.e. such that (h" p )dyn 



-l,a/3 



+ 0(r ). This quantity is a functional of the time varying moments (i.e. having 



^ 2) evaluated at retarded time u — t — r/c, and explicitly reads 



AYi 



£>2 
^>2 



c^+ 2 £! ' 

W1 iL-l\ U ) 



e^2 



2£ 



nL ~ 2 M (e) (u) 



(3.1a) 
(3.1b) 
(3.1c) 



All the logarithms in the quadratic metric h 2 will be generated from the leading 1/r 2 piece in 

~ 2 Qf ( 

4M (2) i k° 



the quadratic source, defined by TV, 



a/3 



from the quantity (3.1) as Q 2 



'z 



+ 



n, u 



+ 0(r ). The coefficient is computed 
a, where the first term will generate the 



tails, and the second term is associated with the stress-energy of gravitational waves, with 

1(1)3^(1)3 _ 1(1)^(1). 



k a = (l,n) the Minkowskian outgoing null vector, and a = z^"^ 1 zi^ — ^> z\^> z\ v . 



Now, as shown in Appendix [AJ the logarithms produced by the second term oc k a k@ are pure 
gauge, so only the first term dealing with tails is responsible for the near-zone logarithms. 
Hence the part of the NZ expansion of the quadratic metric h 2 containing those logarithms 
is given by 



5fi. 



a/3 



+00 



fp y 



k=0 



_d_ 



2k 



A 



-k-l 



B AM (2) 



z{ aP {n,u) 



r 2 c 4 



(3.2) 



We substitute the explicit expression (3.1) into (3.2), expand the retardation u = t — r/c in 
the source term when r — > 0, and integrate using Eqs. (2.9)-(2.10). Then we look for the 



poles oc 1/5 and after applying the finite part get the logarithms. Some general formulas 
for obtaining the logarithms directly from the unexpanded source are relegated to Appendix 
[Aj We readily recover that the dominant logarithms arise at 4PN order. We limit our 
computation to the leading order 4PN and next-to-leading order 5PN logarithms, and find 
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r (6) 
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hi 



(6) 
be 
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(3.3a) 
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+ 



1 

35 



,13 



(3.3b) 
(3.3c) 



The mass- type quadrupole moment M» 3 -, mass octupole moment Mjjfc and current 



quadrupole Sij in Eqs. (3.3) are functions of coordinate time t. The indicated PN remainders 
Ok 



refer only to the logarithmic terms. 



We now want to iterate the expressions (3.3) at higher non-linear order in order to get 



the complete family of logarithms generated by that "seed" . To do that it is very convenient 
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to perform first a change of gauge. Starting from (3.3), which is defined in some harmonic 
gauge, we pose k%f = UZf + 2d^ a ^ — V a ^^n^2 with gauge vector 
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M 
M 



2 x a6 M (5) 1^ a6 M (7) 

3 ab 21c 2 ab 



„abc 
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(3.4a) 



(5) 



In 



(3.4b) 



This gauge transformation will have the effect of moving many 4PN logarithmic terms into 
the 00 component of the (ordinary covariant) metric. As a result the implementation of the 
non-linear iteration in that new gauge will be especially simple. Since our aim is to compute 
the gauge invariant quantity (1.2) we can work in any convenient gauge. Our chosen gauge 
is very similar to the generalization of the Burke-Thorne gauge introduced in [32J to deal 
with higher-order (2.5PN and 3.5PN) radiation-reaction effects. We obtain 
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(3.5a) 
(3.5b) 

(3.5c) 



In this gauge the iteration at cubic non-linear order is very simple. To control all the 
5PN logarithmic terms at cubic order n = 3 we need only to solve the Poisson equation 
A[5F 3 + <5F 3 ] = -20,7? dj5F 2 + C(c- 14 ), where h\ denotes the NZ expansion of the 
linearized metric (2.1a), and we can use for 5k^ the leading 4PN approximation given by 
the first term in (3.5a). Posing h x = —4U/c 2 + (9(c~ 4 ), the latter equation is integrated as 6 



(-7-OO r -j-i 

okn + ok? 



•^ZJx^M^h, 

5 C 12 UX 1V1 ab 111 



o 



(3.6) 



with the explicit expression 



+00 



£=0 



-M L d L 



(3.7) 



We readily check that the quartic and higher non-linear iterations (n ^ 4) are not needed for 
controlling the 4PN and 5PN logarithmic terms (cf. the discussion at the end of Sec. II B). 



6 Actually the integration yields in addition to the near-zone 5PN logarithm (3.6) the extra far-zone 5PN 
logarithmic contribution 
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-12 1V1 ab 111 
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We argued on general grounds in Sec. II C that FZ logarithms do not have to be considered for the present 
computation, so we drop this term out in the following. 
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B. Internal near-zone post-Newtonian metric 



The metric we computed so far is in the form of a multipolar expansion valid in the 
exterior of an isolated source. We now want to deduce from it the metric inside the matter 
source. First of all, since the expressions (3.5) are regular at the origin r — > 0, we find 



using a matching argument that they are necessarily also valid inside the matter source. 



On the other hand it is clear that the expression (3.6) will also be valid inside the source 



provided that we match the multipole expansion U given by (3.7) with the actual Newtonian 



potential of the source. From the known Newtonian limit of the multipole moments Ml = 
f d 3 ii i p(x, t) + 0(c~ 2 ), where p is the Newtonian source density in the source, we get 
GU = U + 0{c~ 2 ) where 



U = G 



dV 



7P( x ' f ) 



(3.8) 



From the latter arguments we therefore obtain the piece of the inner metric of any isolated 
source (coming back to the usual covariant metric g a p) that depends logarithmically on the 
distance r to the source's center at 4PN and 5PN orders as 



(S'tfoo 



G 2 M 
~io 



G 2 M 

c 11 
G 2 M 



,10 



5 

16 
21 



2U 

1? 



x ab M<$ + 



35c 2 



64 

45^ 



Z ioh M% - ~-e iah x M S { 



IIC £j(6) 

be 



g x m ab °v 



111 



o 



In 



,12 



A 



+ o 



189c 2 
1 



x abc M^ c 



In 



A 



+ o 



13 



1 

(3.9a) 
(3.9b) 

(3.9c) 



where U is the Newtonian potential (3.8) valid all over the source 



However we now discuss other pieces of the inner metric whose near-zone expansion 
does not explicitly depend on the logarithms of r but which involve new inner potentials 
integrating over a logarithmically modified source density. The first of these pieces comes 



from the fact that the 4PN modification of the metric given by the first term in (3.9a) implies 



a modification of the stress-energy tensor of the matter fluid at the 5PN order; in particular 
the fluid's source density, say a = T 00 /c 2 , gets modified by the amount 



5a 
P 



4 G 2 M 



5 c 



10 



-x 



ab M^ In 



+ o 



(3.10) 



On the other hand the 4PN term of the metric will induce a 4PN change in the acceleration 
of the fluid motion given by 



5a l 



8G 2 M 



x a M { 3 In 



+ o 
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(3.11) 



When computing the inner metric at the 1PN order we have to take into account the 



retardation due to the propagation of gravity, using say □ —A + c 



~ 2 d 2 A~ 



+ 0(c" 4 ). 



The time derivatives at 1PN order will yield an acceleration and the modification of the 
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acceleration (3.11) will give a contribution at 5PN order. We find that the sum of the two 
effects gives the following extra contribution to the inner metric at 5PN order 



o goo 
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(3.12) 



This 5PN contribution is present only in the 00 component of the metric. 7 The complete 
logarithmic contributions we shall consider in this paper are thus given by 

Sg a /3 = S'gap + S"g a/ 3 . (3.13) 

These contributions exhaust the possibilities of having 4PN and 5PN near-zone logarithmic 
terms in the gauge invariant observable quantity (1.2). 



C. Application to compact binary systems 



Let us now apply the previous results to the specific problem of a system of two point 
particles. The Newtonian mass density in that case is p = ^ a ra a 5(x — Ya) where 5 is 
the Dirac delta function. The trajectory of the a-th particle (a=l,2) is denoted y a (t); the 
ordinary coordinate velocity will be v a = dy a /di. The two masses m a have sum m = 
rrii + m 2, reduced mass pt = rriim^/ 'm and symmetric mass ratio v = pi/m. The Newtonian 
potential of the system reduces to 



U 



Gmi Grrio 
+ 



n 



r 2 



(3.14) 



where r a = |x — y ffl | is the distance from particle a. The regularized value of that potential 
at the location of particle 1 is simply 



_ Gm 2 
u i — ) 



(3.15) 



where r 12 = |yi — y2 1 - Similarly we evaluate the logarithmic contributions at the location of 



particle 1. Concerning the first piece (3.9) we find (no longer mentioning the PN remainder) 
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(3.16a) 
(3.16b) 



7 Interestingly, it was found in Ref. [33] (following [33]) that a similar looking contribution must also be 
taken into account when computing the higher-order (3.5PN) radiation-reaction force for compact binary 
systems from a near-zone radiation-reaction formalism defined in [32]. Actually the 2.5PN+3.5PN near- 
zone radiation-reaction formalism [32j (see in particular Eqs. (2.16) there) is quite similar to the present 
4PN+5PN near-zone conservative logarithm formalism. 
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(3.16c) 



which involves the logarithm ln(r/A) evaluated on the particle 1, i.e. ln(|yi|/A). As for the 
second piece (3.12) we compute the Poisson integral using p = J2 a m a^(' x ~ Ya) an d perform 
a regularization on the particle 1 to obtain 

8G 2 M 
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(3.17) 



which is proportional to the logarithm ln(|y 2 |/A) associated with the other particle. These 
results are valid in a general frame. However we shall later specify the origin of the coordinate 
system to be the center of mass of the binary system. In that case we have ln(|y a |/A) = 
In(ri2/A) + ln(/i/m a ) + 0(c~ 2 ), where the PN remainder does not involve any logarithmic 
term, and the logarithm of the mass ratio is a constant, and is therefore clearly irrelevant 
to our search of logarithmic terms; so In(ri2/A) is in fact the only relevant logarithm and 
we shall now systematically replace all ln(|y a |/A)'s by ln(ri 2 /A). Finally we end up with the 
following contributions of the 4PN and 5PN logarithms to the near-zone metric evaluated 
at the location of particle 1 in our chosen gauge, 
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(3.18a) 
(3.18b) 
(3.18c) 



Note that this result is complete but not fully explicit because we have still to replace all 
the multipole moments and Sl by their expressions valid for point mass binary systems. 
In particular the quadrupole mass moment should be given with 1PN relative precision 
as (1 2 means adding the same terms for particle 2) 
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(3.19) 

and its time derivatives should consistently use the 1PN equations of motion. Besides 
we also need the constant mass monopole or total mass M at 1PN order, namely 



M = 771! 



1 + 



1 



Gm 2 \ 



1 -B-2. 



(3.20) 

(3.21a) 
(3.21b) 

However in applications it is often better to postpone the (messy) replacements of the mul- 
tipole moments by their explicit values (3.19 )— (3.21 ) and to use more compact expressions 
such as (3.18). 



All the other moments are only required at the Newtonian accuracy, and read 

M L = mi y[ + 1 <->• 2 , 
S L = mi e ab{il yi~ l)a v\ + 1 «+ 2 . 
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IV. LOGARITHMS IN THE EQUATIONS OF MOTION AND ENERGY 



General orbits 



With the 4PN and 5PN logarithmic contributions in the near-zone metric (3.18) we 



now derive the corresponding terms in the acceleration of point particle binary systems. 
The computation is straightforward from the geodesic equation. A subtle point is that we 
must take into account the coupling between the 1PN terms in the metric and the 4PN 
logarithm to produce new 5PN logarithms. On the other hand one must be careful about 
the replacement of accelerations in 1PN terms by the 4PN acceleration to also produce 5PN 
logarithms. The final result, valid for generic (non-circular) orbits in an arbitrary frame, is 
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(4.1) 



where the multipole moments are given by (3.19)-(3.21). 



An important check of this result is that the acceleration should be purely conservative, 
by which we mean that there should exist some corresponding contributions at the 4PN and 
5PN orders in the conserved energy, angular momentum, linear momentum and center-of- 
mass position of the binary system. Let us see how this works in the case of the energy. 
The modification at 4PN and 5PN of the energy, say SE, should be such that d5E/dt 
exactly balances the replacement of accelerations by ( |4.1[ ) in the time derivative of the 
known expression of the energy up to 1PN order (say ^ipn)- This requirement yields 
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Plugging (4.1) into (4.2) and using the expressions of the multipole moments (3.19)— (3.21 ), 



we indeed find that the right-hand-side of (4.2) takes the form of a total time derivative, 



and we are thus able to infer the contribution to the energy, 
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In addition to the standard multipole moments ( 3.19[ )-( p.21[ ), we have also introduced the 
supplementary moments (needed only at Newtonian accuracy) 
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By the same method we have also computed the modification of another integral of the 
motion, namely the center-of-mass position G l . The result is 
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and will be useful when restricting our general result for the redshift observable — valid for 
a generic orbit in an arbitrary frame — to circular orbits described in the center-of-mass 
frame. 



B. Circular orbits 



Let us now focus our attention on the case of circular orbits. We look for the 4PN and 



5PN logarithms 5a\ 2 in the relative acceleration a\ 2 



i\ — a 2 of the particles for circular 
orbits. The first contribution to 5a\ 2 will evidently come from the difference Sa\ — 5a 2 . We 
insert the center-of-mass relations y l a = ^[y^, v i2], expressing the individual positions in 
terms of the relative position y\ 2 — y\ — y\ and relative velocity v\ 2 = v \ — v 2 = dy l 12 /dt. At 
1PN order and for circular orbits these expressions simply reduce to the Newtonian relations 
y{ = X 2 y\ 2 and y l 2 = —Xiy{ 2 , where X a = m a /m. All the multipole moments and their 
time derivatives are replaced by their expressions for circular orbits given in terms of y\ 2 and 
v\ 2 (and masses). However there is also another contribution which comes from the known 
relative acceleration at 1PN order (say a^ PN ) when reduced to circular orbits. As usual we 
perform an iterative computation: knowing first 5a\ 2 at 4PN order we use the result to find 
the next order 5PN correction. In this computation we use the fact that the center-of-mass 
relations y l a = ^[y^, v 12 ] are not modified by logarithmic terms before the 5PN order. This 
is checked using the modification of the integral of the center-of-mass 5G % given in (4.5) (see 



also the result (5.1) for circular orbits below, which is clearly a 5PN effect). Finally the 
modification of the acceleration is found to be of the form 
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where the total change in the orbital frequency (squared) for circular orbits due to 4PN and 
5PN logarithms reads 
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The orbital separation is r\% = \yu\, and we have introduced the convenient post-Newtonian 
parameter (where m = mi + m^) 

Gm 
r 12 c 2 



From (4.7 ) we have the relation between the orbital frequency and the parameter 7. Inverting 



this relation we obtain 7 as a function of the orbital frequency or, rather, of the parameter 
x defined by 

'° mn ^ . (4.9) 
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We find that the 4PN and 5PN logarithms in 7 as a function of x are 
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We have taken into account in (4.7) and (4.10) the important fact that the length scale 
A = cP is related to the period P = 2tt/Q, and hence contributes to the logarithm. As 
already mentioned, using Kepler's third law we have = so that ln(-^) = |ln7 
plus an irrelevant constant. Post-Newtonian corrections to Kepler's law do not change the 
argument, which applies with x as well as with 7. Recall that A = 2irc/fl was introduced in 
the problem when we assumed the existence of the helical Killing vector K a d a = d t + Q 
to describe exactly circular orbits. Then this scale entered explicitly into the propagator we 
used to integrate the field [see (2.7) or (3.2)], and it is thus no surprise that it contributes to 
the final result. Of course we could have chosen any other scale proportional to A without 
changing the result which concerns only the logarithmic dependence. 



To be clearer about formulas such as (4.7) and (4.10) we would need to give the more 
complete formulas including also the known contributions up to 3PN order. However we 
must be careful since these formulas depend on the gauge. Thus 5Q 2 and £7 are to be added 
to the 3PN expressions given by Eqs. (188) and (191) in j3] when working in Hadamard 
regularization gauge, or by Eqs. (B6) and (B7) in Paper I when working in dimensional 



regularization gauge. Also the 4PN and 5PN terms computed in (4.7) and (4.10) themselves 



depend on the choice of gauge at the 4PN and 5PN orders (see Sec. III). 

It is much better to turn to gauge invariant quantities. The most obvious one is the 
conserved energy E for circular orbits as a function of the orbital frequency Q. As for the 
previous computation of the acceleration we have two contributions, one coming directly 



from the general-orbit modification of the energy given by (4.3), and one coming from the 
circular-orbit reduction of the 1PN energy Eipn- We first express the results entirely in 



terms of the parameter 7 using (4.7) and then replace all the 7's by functions of the x's 
using (4.10). The result for the 4PN and 5PN logarithms is (where /x = mxm^/ra) 



5E = ~2^ c x 
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Beware that SE here has not the same meaning as in (4.3) because of the additional terms 



coming from the circular-orbit reduction of the 1PN energy £ifn- 

Since the energy as a function of x is a gauge invariant relation, let us also provide 
the complete result including all the known terms up to 3PN order, and also the 4PN 
and 5PN terms in the test-mass limit for one of the particles known from the exact result 
linwo E/(fic 2 ) = (1 - 2x)(l - 3a;)- 1 / 2 - 1. We have 
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(4.12) 



Here eAW) and e^{y) denote some unknown 4PN and 5PN coefficients which are some poly- 
nomials in the symmetric mass ratio v — this can be proved from the fact that the energy 
function for general orbits (i.e. before restriction to circular orbits) must be a polynomial in 
the two separate masses m\ and m2. This 5PN accurate formula could be used to compute 
the location of the innermost circular orbit (ICO) in the comparable mass regime, which also 
coincides with the innermost stable circular orbit (ISCO) in the extreme mass ratio regime. 
The shift of the Schwarzschild ISCO due to the conservative part of the self-force has been 
recently computed |35j. A high-order PN comparison with this result would be interesting, 
but requires at least the evaluation of the coefficients e 4 (z/) and e 5 (z/) in the extreme mass 
ratio regime, i.e. the knowledge of e 4 (0) and e 5 (0). 



V. THE GAUGE INVARIANT REDSHIFT OBSERVABLE 



We are now ready to implement our computation of the gauge invariant redshift ob- 
servable (1.2). We replace the 4PN and 5PN logarithmic terms in the metric coefficients 



evaluated on the particle (3.18) into (1.2). We are careful at including also the metric up 



to 1PN order because of the coupling between 1PN and 4PN orders which produce 5PN 
terms. The result is valid for any orbit in a general frame. Next we go to the frame of the 
center-of-mass defined by G l = 0, where G l is the conserved integral of the center of mass. 
We have found the 4PN and 5PN logarithms in G l in Eq. (4.5), and from this we compute 



the displacement of the center of mass for circular orbits. As already used in Sec. |IV| we find 
that the first logarithmic terms in the center-of-mass integral for circular orbits arise only 
at 5PN order. We obtain 



5G % 
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A 7 5 ln 7? /i 2 , 



(5.1) 



where A = (m2 — m\)jm = yl — 4v is the relative mass difference. The correction to the 
individual center-of-mass positions will thus be given by 5y l a = —5G l /m for a=l,2 (see e.g. 
the Appendix B in Paper I), and similarly 5v l a = —5G l /m for the individual center-of-mass 
velocities. We already notice that because of the factor v 2 in (5.1) this correction will not 



influence the SF limit. Next we reduce the latter expression to circular orbits, replacing all 
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orbital frequencies by their expressions in terms of 7, and then replacing all 7's by their 



expressions in terms of x. The formulas (4.7) and (4.10) for the 4PN and 5PN logarithms 



play of course the crucial role. Finally we end up with the full correction due to the 4PN 
and 5PN logarithmic terms for circular orbits in our redshift observable u T as (removing 
now the index 1 indicating the smaller mass) 
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(5.2) 

This correction is valid for any mass ratio q = m\jmi and is to be added to the 3PN 
expression for u T obtained in Eq. (4.10) of Paper I. Being proportional to the symmetric 
mass ratio u, the correction (5.2) vanishes in the test-mass limit, which is to be expected 



since the Schwarzschild result for u T (Q) does not involve any logarithm. 

We now investigate the small mass ratio regime q <C 1. As in Paper I we introduce a 
convenient PN parameter appropriate to the small mass limit of particle 1: 
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We immediately obtain, up to say the quadratic order in q, and keeping only the relevant 
logarithmic terms, 
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Our complete redshift observable, expanded through post-self-force order, is of the type 
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Schw 



(5.5) 



where the Schwarzschild result is known in closed form as u^ chw = (1 — 3y) ~ 1 ^ 2 : . Adding back 
the 3PN results of Paper I (see Eq. (5.5) there), we thus find that the self-force contribution 
is given by 8 
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(5.6) 



The expansion (5.6) was determined up to 2PN order cx y 3 in [2] based on the Hadamard- 
regularized 2PN metric given in [3]. The result at 3PN order oc y A was obtained in Paper 
I using the powerful dimensional regularization (as opposed to Hadamard's regularization 
which found its limits at that order). By contrast our analytic determination of the logarith- 
mic terms at 4PN and 5PN orders depends only marginally on the regularization scheme. 

The coefficients 0-4 and a$ denote some unknown purely numerical numbers which would 
be very difficult to compute with PN methods, and should depend crucially on having 
a consistent regularization scheme. By comparing the expansion (5.6) with our accurate 



For clarity we add the Landau o symbol for remainders which takes the standard meaning. 
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numerical SF data for «g F , we shall be able to measure these coefficients with at least 8 
significant digits for the 4PN coefficient 0:4, and 5 significant digits for the 5PN coefficient 
a 5 . These results, as well as the estimation of even higher-order PN coefficients, will be 
detailed in Sec. IVT1 

Similarly, adding up the results of Paper I for the post-self-force term, we get 

rp 2 97 o /725 41 2 \ 4 

+ e 4 y 5 + (e 5 + \ny)y e + o(y 6 ) . (5.7) 



35 

Note that there is no logarithm at 4PN order in the post-self-force term, as is also seen from 



Eq. (5.4); the next 4PN logarithm would arise at cubic order q , i.e. at the post-post-SF 



level. The coefficients €4 and €5 in (5.7) are unknown, and unfortunately they are expected 
to be extremely difficult to obtain, not only analytically in the standard PN theory, but also 
numerically as they require a second-order perturbation SF scheme. 

VI. NUMERICAL EVALUATION OF POST-NEWTONIAN COEFFICIENTS 

In the self-force limit, the SF effect u^ F on the redshift observable u T is related to the 
regularized metric perturbation h^g at the location of the particle through 

«sF = ^(l-3y)- 1/2 ^^, (6.1) 

where u a is the background four-velocity of the particle. Beware that here stands in fact 
for the perturbation per unit mass ratio, denoted h^/q in Paper I (cf. Eq. (2.11) there). In 
SF analysis, the combination u a u l3 h^ l3 arises more naturally than Ug F ; this is the quantity 
we shall be interested in fitting in this Section. However our final results in Table [V] will 
include the corresponding values of the coefficients for the redshift variable «g F . We refer to 
Sec. II of Paper I for a discussion of the computation of the regularized metric perturbation 
h^Ri an d the invariant properties of the combination W^u^h^^ with respect to the choice of 
perturbative gauge. In this Section we often use r = 1/y, a gauge invariant measure of the 



orbital radius scaled by the black hole mass 1x12 [see Eq. (5.3) . 

Our earlier numerical work, partially reported in [2] and in Paper I, provided values of 
the function u a u l3 h^ j3 (r) which cover a range in r from 4 to 750. Following a procedure 
described in we have used Monte Carlo analysis to estimate the accuracy of our values 
for u a u^h^g. As was reported in Paper I, this gives us confidence in these base numbers to 
better than one part in 10 13 . We denote a standard error a representing the numerical error 
in u a u^h^p by 

a ~ \u a u fi h%\ x E x 1(T 13 , (6.2) 

where E ~ 1 is being used as a placeholder to identify our estimate of the errors in our 
numerical results. 

A. Overview 

A common task in physics is creating a functional model for a set of data. In our problem 
we have a set of iV data points /j and associated uncertainties a iy with each pair evaluated 
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at an abscissa r^. We wish to represent this data as some model function f(r) which consists 
of a linear sum of M basis functions Fj(r) such that 

M 

The numerical goal is to determine the M coefficients Cj which yield the best fit in a least 
squares sense over the range of data. That is, the Cj are to be chosen such that 



N 



* 2 = £ 



i=i 



M 

, el' A r, i 

(6.4) 



is a minimum under small changes in the Cj. For our application we choose the basis functions 
Fj(r) to be a set of terms which are typical in PN expansions, such as r -1 , r~ 2 , . . . , and 
also terms such as r _5 ln(r). 

Our analysis depends heavily upon Ref . [37] ; we use both the methods and the computer 
code for solving systems of linear algebraic equations as described therein. While we do 
employ standard, least squares methods for solving a system of linear equations to determine 
the Cj, we also recognize that a solution to this extremum problem is not guaranteed to 
provide an accurate representation of the data (fj, /j, <7i). The quality of the numerical fit 



is measured by y 2 as defined in Eq. (6.4). If the model of the data is a good one, then 
the x 2 statistic itself has an expectation value of the number of degrees of freedom in the 
problem, N — M, with an uncertainty (standard deviation) of >J2(N — M). In particular, a 
large residual x 2 would correspond to under-fitting the data whereas a x 2 that is too small 
corresponds to over-fitting the data, which amounts to fitting randomness in the residuals. 

The numerical evaluation of the fitting coefficient Cj includes a determination of its un- 
certainty Ej which depends upon i) the actual values of T\ in use, ii) all of the <7j, and iii) the 
set of basis functions Fj(r). In fact, the estimate of the Ej depends solely upon the design 
matrix 

A l3 = ^ , (6.5) 

and not at all on the data (or residuals) being fitted. However, the estimates of the Ej are 
only valid if the data are well represented by the set of basis functions. For emphasis: the 
Ej depend upon Fjfa) and upon Oi but are completely independent of the fi. Only if the fit 
is considered to be good, could the Ej give any kind of realistic estimate for the uncertainty 
in the coefficients Cj. If the fit is not of high quality (unacceptable x 2 ), then the Ej bear no 
useful information [37]. We will come back to this point in the discussion below. 

A further remark concerning the meaning of the Ej is appropriate. Fitting the data as 
described to determine the coefficients is a standard, well defined statistical procedure. If 
we were to change the integration routine used to generate the u a vP h^Jr^ ; which are the 
set of input data values fi, with the restriction that we maintained the same numerical 
accuracy then the fi would each change in a random way governed by <jj. If the coefficients 
Cj were then determined for this second data set, the statistical analysis ensures that the 
Ej associated with this second data set would be identical to those of the first set and the 
newly determined estimates of the Cj would differ from the initial ones in a statistical fashion 
governed by the Ej. Recall that the Ej depend upon the choice of the r i; upon the <7j for 
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the individual data points and upon the set of basis functions Fj(r). The Ej are completely 
independent of the data values /j. 

Now we consider two other possible changes. If we add an extra data point, or if we add 
another basis function not orthogonal to the others (this would be typical over a finite set of 
data points, unless we carefully engineered otherwise) the design matrix changes accordingly, 
all estimated coefficients Cj change accordingly, and the estimated Ej change in ways which 
are not easily related to the previous results. In particular, if we add an additional basis 
function Fm+i to the previous set, so there is now one more coefficient Cm+i to be fit, and 
we compare the first M values of the new Cj to their earlier values, their differences need 
not be closely related to either the first or second set of Ej. Thus, a change in the design 
matrix of the problem leads to an inability to make any intuitive prediction about what to 
expect for the new Cj, and there is no reason to expect that the differences of the Cj respect 
the values of the E^ for these two different statistical problems. 

We also should remark that the task of determining coefficients in the 1/r characterization 
of our numerical data is almost incompatible with the task of determining an asymptotic 
expansion of u a Wh^g from an analytic analysis. Analytically, the strict r — > +00 limit 
is always technically possible, whereas numerically, not only is that limit never attainable, 
but we must always contend with function evaluations at just a finite number of discrete 
points, obtained within a finite range of the independent variable, and computed with finite 
numerical precision. Nevertheless, this is what we intend to do. 

In practice, the numerical problem is even more constrained. At large r, even though the 
data may still be computable there, the higher order terms for which we are interested in 
evaluating PN coefficients rapidly descend below the error level of our numerical data. This 
is clearly evident in Fig. [T] below. For small r, the introduction of so many PN coefficients is 
necessary that it becomes extremely difficult to characterize our numerical data accurately. 
Thus, in practice, we find ourselves actually working with less than the full range of our 
available data. At large r we could effectively drop points because they contribute so little 
to any fit we consider. At the other extreme, the advantage of adding more points in going 
to smaller r is rapidly outweighed by the increased uncertainty in every fitted coefficient. 
This results from the need to add more basis functions in an attempt to fit the data at small 
r. Further details will become evident in Sec. IVI Dl below. 



B. Framework for evaluating PN coefficients numerically 

In a generic fashion we describe an expansion of u a u l3 h^ /3 in terms of PN coefficients a 



and bj with 



z-.a-.l 



where ao is the Newtonian term, a\ is the 1PN term and so on. Similarly, for use in 
applications involving u T we also introduce the coefficients aj and (3j in the expansion of 
the SF contribution 

/ -^^--InrV^i- 

These series allow for the possibility of logarithmic terms, which are known not to start 
before the 4PN order. We also concluded in Sec.[n]that (lnr) 2 terms cannot arise before the 
5.5PN order. Since we are computing a conservative effect, possible time-odd logarithmic 
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squared contributions at the 5.5PN or 6.5PN orders do not contribute. But there is still the 
possibility for a conservative 7PN (lnr) 2 effect, probably originating from a tail modification 
of the dissipative 5.5PN (lnr) 2 term. However, we shall not permit for such a small effect in 



our fits. As discussed below in Sec. VI D we already have problems distinguishing the 7PN 
linear lnr term from the 7PN non-logarithmic contribution. 

The analytically determined values of the coefficients ao, a\, a 2 , and ao, ai, a 2 , «3 
computed in Ref. [2] and Paper I are reported in Table [l| together with the new results 
64 = -*f , h = fnf and & = -f , & = fi of the present work. 
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TABLE I: The analytically determined PN coefficients for u a u^h^a (left) and Ug F (right). 



C. Verifying analytically determined PN coefficients 

In this Section we investigate the use of our data for u a u l3 h^ l3 and the fitting procedures 



we have described above (and expanded upon in the beginning of Sec. VI D). We will begin 
by fitting for enough of the other PN coefficients to be able to verify numerically the various 
coefficients 03, 64 and 65 now known from PN analysis. We choose a starting point for the 
inner boundary of the range, and each range continues out to r = 700. The results of a series 
of fits are displayed in Tables |TT] and III First we remark that bringing the outer boundary 



inward as far as to 300 has very little effect on the outcome of any of these fits, except that 
the x 2 statistic decreases as expected with the number of degrees of freedom. 

As a first step in this Section, we will complete the task we began in Paper I, namely, the 
numerical determination of the coefficient a 3 (and a 3 ), this time taking fully into account 
the known logarithmic terms at 4PN and 5PN order. For illustrative purposes only, these 
results are given in Table [IT] We were able to obtain a fit with six undetermined parameters, 
and could include data from r = 700 down to r = 35. Note that, with the inclusion of the 
64 and 65 coefficients, the precision of our tabulated value for 03 has increased by more than 
four orders of magnitude from Paper I, although our accuracy is still no better than about 
2E. Such a discrepancy is not uncommon. The uncertainty, E, reflects only how well the 
data in the given, finite range can be represented by a combination of the basis functions. It 
is not a measure of the quality of a coefficient when considered as a PN expansion parameter, 
which necessarily involves an r — > +00 limiting process. 

Our next step is to include the known value for 03 and to use our numerical data to 
estimate values for the 64 and b 5 coefficients. Our best quality numerical result was obtained 
with five fitted parameters, over a range from r = 700 down to only r = 65, and is given 



in the first row of Table III Notice that while our 6 4 is determined relatively precisely, it 



24 



rr 

coeft. 


value 


a 3 


-32.5008069(7) 


(24 


— 1z1.oUZ04^oU ) 


a 5 


-42.99(5) 
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TABLE II: The results of a numerical fit for a set of coefficients which includes the analytically 
known 03. Thus this is not the best-fit of our data possible, but it allows for a comparison with 
Table |VJ The uncertainty in the last digit or two is in parentheses. The range runs from r = 35 to 
r = 700, with 266 data points and a respectable % 2 of 264. 



has only about 6E accuracy. The higher order coefficient 65 is more difficult to obtain and, 
at this point, it is very poorly determined. It corresponds to a term which falls off rapidly 
with increasing r and is significant over a relatively small inner part of the fitted range. 

We can of course use the known value of 64 in order to improve the accuracy for 65. If 
we do this without adding another parameter to fit, we immediately get a fit of very poor 
quality, since we have moved 64 far from its best-fit value; as shown in the second row of 
Table [TTll we must move the inner boundary out to r = 85 to re-establish a good fit. 
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TABLE III: The numerically determined PN coefficients for u a u^h^a. Each row represents a 
different fit. The first three columns give the starting point r m ; n at the inner boundary of the 
fitting range, the degrees of freedom of the fit, N — M, and the x 2 statistic for the chosen fit. If a 
value for a coefficient is not shown, then either that parameter was not included in that particular 
fit (far right) or its analytically known value was used (e.g., 64). The formal uncertainty of a 
coefficient in the last digit or two is in parentheses. The outer boundary is 700 in each case. 



The inclusion of basis functions for the higher order coefficients, b$ and 07, as shown in 
the third and fourth rows, respectively, allows the inner boundary for the fit to move to 
smaller r where the higher PN terms are more important. The third row of the table shows 
that adding another parameter allows us to move the inner boundary to r = 65, while the 
final row shows that we can now add one further fitted parameter, and obtain a good quality 
fit by pushing the inner boundary to r = 40. Only in this row is the 65 parameter close to 
its known value, but it is still off by around 4.5S (see Table IV below). Moreover, we have 
reached a limit for treating our data in this way, adding further parameters and inner points 
does not result in any higher quality fit. 

By now we have presented enough to show that we have data which allows high precision, 
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with an accuracy that we now have some experience in relating to the computed error 
estimates. This experience will be valuable when we come to discuss further results in the 



next Section. For convenience, we summarize the relevant information further, in Table IV 



referring just to our estimates of known PN parameters, and relating our error estimates to 
the observed accuracy. 
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TABLE IV: Comparing the analytically known PN coefficients (column 5) with their numerically 
determined counterparts (column 3), and comparing the numerically determined error estimates 
(column 3) with the apparent accuracy (column 4). The source of the data is given in column 1. 



D. Determining higher order PN terms numerically 

In this Section we turn our attention to using our numerical SF data and fitting procedures 
to obtain as many as possible unknown PN coefficients, by making maximum use of the 
coefficients which are already known. We find that in our best fit analysis we can use a set 
of five basis functions corresponding to the unknown coefficients 04, 05, a 6 , 6 6 and 07. 
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TABLE V: The numerically determined values of higher-order PN coefficients for u a u^h^^ (left) 
and for -Ug F (right). The uncertainty in the last digit or two is in parentheses. The range runs 
from r = 40 to r = 700, with 261 data points being fit. The x 2 statistic is 259. We believe that a 
contribution from a 67 term piggybacks on the aj coefficient. Both terms fall off rapidly and have 
influence over the fit only at small r. And the radial dependence of these two terms only differ by 
a factor of lnr [or possibly (lnr) 2 ] which changes slowly over their limited range of significance. 



In Table [V] we describe the numerical fit of our data over a range in r from 40 to 700. The 
X 2 statistic is 259 and slightly larger than the degrees of freedom, 256, which denotes a good 
fit. Further, we expect that a good fit would be insensitive to changes in the boundaries 
of the range of data being fit, and we find, indeed, that if the outer boundary of the range 
decreases to 300 then essentially none of the data in the Table changes, except for x 2 an d 
the degrees of freedom which decrease in a consistent fashion. Figure [T] shows that in the 
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FIG. 1: The absolute value of the contributions of the numerically determined post-Newtonian 
terms to r 5 u a u^h^ l3 . Here PNL refers to just the logarithm term at the specified order. The 
contribution of 04 is not shown but would be a horizontal line (since the 4PN terms behaves like 
r~ 5 ) at approximately 121.3. The remainder after 04 and all the known coefficients are removed 
from r 5 u a u^h^p is the top (red) continuous line. The lower (black) dotted line labelled "err" shows 
the uncertainty in r 5 u a u /3 h^ /3 , namely 2Er 4 x 10~ 13 . The jagged (green) line labelled "|res|" is the 
absolute remainder after all of the fitted terms have been removed. The figure reveals that, with 
regard to the uncertainty of the calculated u a u^h^ l3 , the choice E ~ 1 was slightly too large. 

outer part of the range u a Wh^g is heavily dominated by only a few lower order terms in the 
PN expansion — those above the lower black double-dashed line in the figure. 

The inner edge of the range is more troublesome. The importance of a given higher 
order PN term decreases rapidly with increasing r. Moving the inner boundary of the range 
outward might move a currently well determined term into insignificance. This could actually 
lead to a smaller x 2 , but it would also lead to an increase in the Ej of every coefficient. 
Moving the inner edge of the range inward might require that an additional higher order 
term be added to the fit. This extra term loses significance quickly with increasing r so 
the new coefficient will be poorly determined and also result in an overall looser fit with an 
increase of Ej for all of the coefficients. If the inner boundary and the set of basis functions 
are chosen properly, then a robust fit is revealed when the parameters being fit are insensitive 
to modest changes in the boundaries of the range. The fit described in Table [V] appears to 
be robust. The parameters in this Table are consistent with all fits with the inner boundary 
of the range varying from 35 to 45 and the outer boundary varying from 300 to 700. 

If an additional term, with coefficient 67, is added to the basis functions then, for identical 
ranges, each of the E 3 - increases by a factor of about ten, and the changes in 04 and 05 are 
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within this uncertainty. The coefficient ae changes sign and be and a 7 change by an amount 
significantly larger than the corresponding E^. And the new coefficient b 7 is quite large. 
In the context of fitting data to a set of basis functions these are recognized symptoms of 
over-fitting and imply that the extra coefficient degrades the fit. 

How should we (and others) interpret the data in Table M? To guide our discussion of 



this very important question, we assemble together into Table VI all the relevant results 



from the earlier fits of Sec. |VI C| which relate to the best prior estimates we have there for 
a>4, 0-5, be and a 7 which we have finally calculated here. As was shown in Table IV and 



is now evident in Table [VI] our numerical accuracy tends to be in the range of 2 — 6£, both 
when comparing the best results for a 4 , a 5 , a 6 , 6 6 and a 7 from Sec. |VI C|with those obtained 



here and, we would suggest, for the purposes of comparing the results of this Section with 
future PN coefficients. 



coeff. 


Tabic 


V 


(best) 


Table 


II 




Table 


III 




04 
a 5 
a 6 

be 
a 7 


-121.30310(10) 

-42.89(2) 

-215(4) 

+680(1) 

-8279(25) 


-121.30254(30) -»• (56) 
-42.99(5) -> (10) 
-228(6) -> (13) 
+677(2) -> (3) 
-8226(27) -»■ (53) 


-121.3052(6) -> (21) 
-47(1) -> (4) 
-359(41) (144) 
+625(15) -»■ (55) 
-7722(162) -> (557) 



TABLE VI: Comparing the "best fit" numerical values and statistical uncertainties of the estimated 



PN coefficients in Table \V\ to other numerical evaluations of these same quantities in Sec. VIC 



E. Summary 

Our best fit can be visualized in Fig. [2J where we plot the self-force effect «g F on the 
redshift variable function of r = y 1 , as well as several truncated PN series up to 

7PN order, based on the analytically determined coefficients summarized in Table |TJ as well 
as our best fit of the higher-order PN coefficients reported in Table |Vj Observe in particular 
the smooth convergence of the successive PN approximations towards the exact SF results. 
Note, though, that there is still a small separation between the 7PN curve and the exact 
data in the very relativistic regime shown at the extreme left of Fig. |2j 

We have found that our data in the limited range of 35 ^ r ^ 700 can be extremely 
well characterized by a fit with five appropriately chosen (basis) functions. That is, the 
coefficients in Table [V] are well determined, with small uncertainties, and small changes in 
the actual details of the fit result in coefficients lying within their error estimates. Fewer 
coefficients would result in a very poor characterization of the same data while more coef- 
ficients result in large uncertainties in the estimated coefficients, which themselves become 
overly sensitive to small changes in specific details (such as the actual choice of points to 
be fitted). In practice, over the data range we finally choose, and with the five coefficients 
we fit for, we end up with exceedingly good results for the estimated coefficients, and with 
residuals which sink to the level of our noise. We have a very high quality fit which is quite 



insensitive to minor details. Nevertheless, as Tables IV and VI hint, error estimates for these 



highest order coefficients should be regarded with an appropriate degree of caution. 
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FIG. 2: The self- force contribution u§ F to u T plotted as a function of the gauge invariant variable y 1 . 



Note that y 1 is an invariant measure of the orbital radius scaled by the black hole mass m 2 [see Eq. (5.3)] 



The "exact" numerical points are taken from Ref. [2]. Here, PN refers to all terms, including logarithms, 
up to the specified order (however recall that we did not include in our fit a log-term at 7PN order). 
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Appendix A: Formulas to compute the PN logarithms 



In Sec. |III| we looked for poles generating near-zone logarithms when integrating the field 
equations at quadratic non-linear order. We used the propagator of the "instantaneous" 
potentials defined by 

+ 00 / r\ \ 2k r> 

fc=0 v 7 



and acting on a source term of the type r~ 2 F(n, u) where u = t — r/c; see Eq. (3.2). We 
consider here a single multipolar piece in the source term, say r~ 2 hiF{u). The function F 
is typically a product of the mass with some time derivatives of multipole moments. We 



recall that the propagator (Al) depends on the length scale A = cP, where P is the period 
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of the source; we thus consider 



F(u) 



(A2) 



In this Appendix we shall provide a general and compact formula giving all the logarithms 
in the near- zone expansion of the solution (A2). The logarithms come from expanding the 
retardation u = t — r/c in the source when r/c — > 0, integrating each of the terms using 



the formulas (2.9)-(2.10), and finally taking the finite part (FP) associated with the poles 
oc B^ 1 . Our compact formula gives the result of all these operations as 



Mi 



In 



d. 



_ r/ c ) - F(- £ - l \t + r/c) 



(A3) 



where _F (_£_1 ' ) denotes the (£+l)-th time anti-derivative of the function F. By 5&l we mean 
the contribution of logarithms in thus all the other terms in besides 5§l admit an 
expansion when r — > in simple powers of r without logarithms. Note that the factor of 



the logarithm in Eq. (A3) is a multipolar antisymmetric homogeneous solution of the wave 



equation which is regular at the origin, when r — > 0. The logarithms in (A3) are thus of 
the NZ type; no FZ logarithms are generated from a source term r~ 2 niF(u). We recall also 
from Sec. [TTJthat the FZ logarithms start to arise at the cubic n = 3 non linear iteration, and 
that they do not contribute to the conservative part of the dynamics of compact binaries. 
The formal near-zone expansion of 5$^ reads 



8®, 



'-Yin 



E 

i=0 



til r 



2i+£ 



2H\(2i + 2£ + c 2i+e 



(A4) 



At the 1PN relative order required for our computation in (3.2)-(3.3), we have 

^,2 



5$, 



(2£ + l)\\c e 



F«\t) 



2c 2 (2£ + 3) 



F^ +2 \t) + 



1 



c> 



In 



A 



(A5) 



The result (|A3j) can be generalized in the following sense that the same type of result will 
hold also for non-STF sources. Namely, if we define 5<&l to be non-STF in L, i.e. having 

— n h ' ' ' n ii i n place of the STF product h L in (A2), then we can easily prove that the 
log-terms are given by (A3) with = d i± ■ ■ ■ d ie in place of the STF product d^. Of course 



all the other terms will be different, but the structure of the log-terms will be the same. 
Then it is trivial to show that the formula applies as well to a product of Minkowskian 
outgoing null vectors k a = (1, n) representing the direction of propagation of gravitational 
waves, and satisfying r\ a pk a k^ = 0. Considering 



OLl'~OL£ 



X~ 



k ■ ■ ■ k 



-F(u) 



(A6) 



where k a = (— l,n), we find indeed that the contribution of logarithms in the near-zone 
expansion of this object is given by 



5$ 



-c 



Ml 



In ( - ) d ai . ..,„ 



Fi-'-Vfr -r/c)- F^- l \t + r/c) 



(A7) 
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We use this result to show that a family of logarithms not considered in Sec. Ill is actually 
pure gauge. We showed there that all the 4PN and 5PN near-zone logarithms come from 



iterating the leading-order 1/r part of the quadratic source, namely Q 2 



a/3 



4M (2) a/3, 

(A *1 



a. 



However we computed only the first term oc 



(2) *P 



which is associated with tails, but we left 



out the second term oc k a W . Now thanks to the structure cx k a hP the logarithms appear 
in the form of a gauge transformation and will never contribute to a gauge invariant result. 
This was already shown at the level of the dominant 4PN log-term in [25] . By expanding a 
on the basis of STF tensors Ul (or rather n^_ 2 ) we need only to prove this for each of the 
individual multipolar pieces in the source which have the structure 



4> 



a/3 
L-2 



T 



-i 



k a k p 



n L -2 



F(u) 



(A8) 



Applying (A7) the logarithms are given by 

\t+i 



-c 



In ( r -) d^d^ ( ~ T/C) ~ F{ ~ £ ~ 1){t + T/C) 



(A9) 



and can readily be put in the form of a gauge transformation with gauge vector 

_ r /c) - F(-«- 1 )(t + r/c) 



L-2 



\nl-)d a d L - 2 



(A10) 



Indeed we have 53>"_ 2 



2d (a ^ L l 2 - v a(3 d^_ 2 modulo some terms which are free of log- 
arithms. Therefore the "seed" logarithms generated in this way at quadratic order can be 
removed by a gauge transformation, and we conclude that the whole family of logarithms 
coming from the iteration at cubic and higher orders can be removed by a non-linear defor- 
mation of the gauge transformation, namely by a coordinate transformation. Thus we do 
not have to consider these logarithms in our computation of a gauge invariant quantit y; on ly 



those coming from the first term cx 



in Q 2 will contribute as computed in Sec. 
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